


% updated version of kymtrack; for very high quality kymographs, this
% actually performs worse than kymtrack_v2_0 since this doesn't have a
% function fitting step. But it is much more robust to noise, movement
% artifacts, and extraneous signals in the frame. 
% 
% Down the line, can add a fitting function to do the same kind of fit as
% kymtrack_v2_0 does, coupled with the added functionality in this
% function. But since doing many linescans per vessel anyway, the averaging
% should take care of any shortcomings
%
% refineposition is either true or false; if true, then additional
% refinement will be done on the position of the edges by fitting them to a
% sigmoid

function [r,rrefine] = kymtrack_v3_4(kymo,refineposition)

%% gather properties:
numframes = size(kymo,2);

%% set some parameters for the function:
doplot = false;
if nargin < 2 || isempty(refineposition)
    refineposition = false;
end
baselinepercent = 7;
buffer = 3;     % number of pixels of buffer to grant the estimated edges (crops signal outside the buffer)
ridgePositionChangePenalty = .75;    %penalty for changing the position of a ridge over time (avoids jumps)
centerCrossingThresh = round(numframes*.1);     %threshold for maximum number of times the edge estimate is allowed to cross the centerline estimate
baselineInten = prctile(kymo(:),baselinepercent);
refinePositionBuffer = [11,4];

%smoothing parameters
kernx = (-8:8)';
sigG = 2;
sigM = 2;        %for mexican hat kernel


%% smooth kymograph IN SPATIAL DIMENSION ONLY:
kymo(isnan(kymo)) = baselineInten;
kernG = exp(-(kernx.^2)/(2*sigG)^2);
kernG = (kernG/sum(kernG(:)));
% generates a mexican hat kernel along the spatial dimension ONLY
kern = 2/(sqrt(3*sigM)*pi^(.25))*(1-(kernx/sigM).^2).*exp(-(kernx.^2)/(2*sigM^2));

% specific kernl to emphasize BOTTOM edge 
kern_btm = kern;
kern_btm(kernx<0) = kern(kernx==0);
kern_btm = kern_btm/sum(kern_btm);

% specific kernl to emphasize TOP edge 
kern_top = kern;
kern_top(kernx>0) = kern(kernx==0);
kern_top = kern_top/sum(kern_top);


smimTop = conv2(kymo,kern_top,'same');
smimTop(smimTop<0) = 0;
%normalize the smoothed kymograph:
smimTop = smimTop-min(smimTop(:));
smimTop = smimTop/max(smimTop(:));
thresh = multithresh(smimTop,2);
segTop = imquantize(smimTop,thresh); 

smimBtm = conv2(kymo,kern_btm,'same');
smimBtm(smimBtm<0) = 0;
%normalize the smoothed kymograph:
smimBtm = smimBtm-min(smimBtm(:));
smimBtm = smimBtm/max(smimBtm(:));
thresh = multithresh(smimBtm,2);
segBtm = imquantize(smimBtm,thresh); 
segIm = min(segTop,segBtm);
temp = segIm > 1;
temp = imfill(temp,'holes');
ind = find(segIm == 1 & temp == 1);
segIm(ind) = 2;

%% estimate centerline and edges of the kymograph signal
temp = min(smimTop,smimBtm);
[centerline,edges] = fitcenterline(temp);
centerline = smooth(centerline,'rlowess');
edges(:,1) = round(smooth(edges(:,1),'rlowess'));
edges(:,2) = round(smooth(edges(:,2),'rlowess'));

%% crop out signal outside the edges of the signal
cropsmim = smimTop;
minimval = prctile(cropsmim(:),10);
for n = 1:numframes
    a = edges(n,1)-buffer;
    b = edges(n,2)+buffer;
    if (a<1)
        a = 1;
    end
    if b > size(kymo,1)
        b = size(kymo,1);
    end
    cropsmim(1:a,n) = cropsmim(a,n);
    cropsmim(b:end,n) = cropsmim(b,n);
    ind = find(segIm(:,n)==1);
    ind = ind(ind < centerline(n));
    cropsmim(ind,n) = minimval;
end
smimTop = cropsmim;

cropsmim = smimBtm;
minimval = prctile(cropsmim(:),10);
for n = 1:numframes
    a = edges(n,1)-buffer;
    b = edges(n,2)+buffer;
    if (a<1)
        a = 1;
    end
    if b > size(kymo,1)
        b = size(kymo,1);
    end
    cropsmim(1:a,n) = cropsmim(a,n);
    cropsmim(b:end,n) = cropsmim(b,n);
    ind = find(segIm(:,n)==1);
    ind = ind(ind > centerline(n));
    cropsmim(ind,n) = minimval;
end
smimBtm = cropsmim; 

if doplot
    figure
    subplot(2,1,1)
    imagesc(kymo);
    subplot(2,1,2)
    imagesc(min(smimTop,smimBtm));
end

%% get edge images of the kymographs:
numridges = 2;
[~, Gy_top] = gradient(smimTop);       %calculate the gradient along the linescan; emphasizing top edge
[~, Gy_btm] = gradient(smimBtm);       %calculate the gradient along the linescan; emphasizing bottom edge

Gy_btm(Gy_btm>0) = 0;
Gy_top = -Gy_top;
Gy_top(Gy_top>0) = 0;
% calculate position of bottom edge
btmRidgeIm = makeridge(Gy_btm);
btmRidges = tfridge(btmRidgeIm,1:size(btmRidgeIm,1),ridgePositionChangePenalty,'NumRidges',numridges);

%temp = centerline + (edges(:,2)-centerline)/2;
temp = centerline;
temp = btmRidges-repmat(temp,1,numridges);
btmRidges(temp<0) = nan;
temp = sum(isnan(btmRidges),1);
if min(temp)>centerCrossingThresh
    btmRidges = max(btmRidges,[],2);
else
    btmRidges(:,temp>centerCrossingThresh) = nan;
    [~,I] = min(btmRidges,[],2);
    ind = mode(I);
    btmRidges = btmRidges(:,ind);
    temp = btmRidges - centerline;
    btmRidges(temp<0) = centerline(temp<0);
end
% fill in empty positions with interpolation:
temp = smooth(btmRidges,'rlowess');
btmRidges(isnan(btmRidges)) = round(temp(isnan(btmRidges)));

% calculate position of top edge
topRidgeIm = makeridge(Gy_top);
topRidges = tfridge(topRidgeIm,1:size(topRidgeIm,1),ridgePositionChangePenalty,'NumRidges',numridges);

temp = centerline;
temp = topRidges-repmat(temp,1,numridges);
topRidges(temp>0) = nan;
temp = sum(isnan(topRidges),1);
if min(temp)>centerCrossingThresh
    topRidges = min(topRidges,[],2);
else
    topRidges(:,temp>centerCrossingThresh) = nan;
    [~,I] = max(topRidges,[],2);
    ind = mode(I);
    topRidges = topRidges(:,ind);
    temp = centerline-topRidges;
    topRidges(temp<0) = centerline(temp<0);
end
% fill in empty positions with interpolation:
temp = smooth(topRidges,'rlowess');
topRidges(isnan(topRidges)) = round(temp(isnan(topRidges)));

%% refine edge positions and return the value:
r = nan(2,numframes);
temp = topRidges;
temp(isnan(temp)) = edges(isnan(temp),1);
r(1,:) = roughEdgeRefine(smimTop,temp,-1);
temp = btmRidges;
temp(isnan(temp)) = edges(isnan(temp),2);
r(2,:) = roughEdgeRefine(smimBtm,temp,1);


if refineposition
    rrefine = nan(size(r));
    for n = 1:numframes
        ind = (r(1,n)-refinePositionBuffer(1)):(r(1,n)+refinePositionBuffer(2));
        ind = ind(ind>0);
        y = kymo(ind,n);
        try
            temp = refineEdgePosition_v1_2(y);      %previous version used refineEdgePosition instead
        catch
            disp('missed a frame')
            continue;
        end
        rrefine(1,n) = ind(1)+(temp(2)-1);
        
        ind = (r(2,n)+refinePositionBuffer(1)):-1:(r(2,n)-refinePositionBuffer(2));
        ind = ind(ind<=size(kymo,1));
        y = kymo(ind,n);
        try
            temp = refineEdgePosition(y);
        catch
            disp('missed a frame')
            continue;
        end
        rrefine(2,n) = ind(1)-(temp(2)-1);
    end
else
    rrefine = [];
end

if doplot
    figure;imagesc(btmRidgeIm);title('bottom ridges')
    figure;imagesc(topRidgeIm);title('top ridges')
    figure
    imagesc(kymo);
    hold on
    plot(topRidges,'r.')
    plot(btmRidges,'g.')
    plot(r(1,:),'r*')
    plot(r(2,:),'g*')
end

%takes the kymograph of the FILLED VESSEL and approximates a center line
%from which you can get the top/bottom components
function [centerline,edges] = fitcenterline(im)
%first frame is garanteed to be centered because that's how the selection
%is done:
intenthresh = .3;
numframes = size(im,2);
centerline = nan(numframes,1);
centerline(1) = round(size(im,1)/2);
edges = nan(numframes,2);
w = exp(-((5:-1:0).^2)./(2*2^2));w = w/sum(w);
for n = 1:numframes
    if n == 1
        est = centerline(n);
    elseif n < (length(w)+1)
        est = centerline(n-1);
    else
        temp = centerline((n-length(w)):(n-1));
        est = round(w*temp);
    end
    top = im(est:-1:1,n);
    btm = im((est+1):end,n);
    top = find(top<(top(1)*intenthresh),1);
    btm = find(btm<(btm(1)*intenthresh),1);
    try
        edges(n,1) = est-top;
        edges(n,2) = est+btm;
    catch
        disp('oh no - centerline finding error!')
        edges(n,1) = edges(n-1,1);
        edges(n,2) = edges(n-1,2);
    end
    
    centerline(n) = round(mean(edges(n,:)));
end
return;


% finds ridges in the image from the eigenvalues of the Hessian
function ridgeim = makeridge(im)

e2thresh = .00007;        %works for gradient; not for imgradientxy
e1thresh = .5;

im = imgaussfilt(im,3);
[gx,gy] = gradient(im);
[gxx,gxy] = gradient(gx);
[~,gyy] = gradient(gy);

e1 = zeros(size(im));
e2 = zeros(size(im));
for n = 1:numel(e1)
    temp = eig([gxx(n),gxy(n);gxy(n),gyy(n)]);
    e1(n) = temp(1);
    e2(n) = temp(2);
end
% want areas where e2 is high and e1 is close to zero
mask = ones(size(e2));
mask(e2<e2thresh) = 0;
mask(e1>(e1thresh*e2)) = 0;
ridgeim = mask.*e2;
return;



function ep = refineEdgePosition_local(im,edgeestimate)
%fits a gaussian to the edge image 
numframes = size(im,2);
ep = edgeestimate;
buff = 10;  % number of pixles on each side to err by
x = 1:size(im,1);x = x';
po = nan(1,3);
opts = optimset('display','off');
for n = 1:numframes
    y = im(:,n);
    y = y-min(y);
    y(1:(edgeestimate(n)-buff)) = 0;
    y((edgeestimate(n)+buff):end) = 0;
    if max(y) == 0
        continue;
    end
    F = @(p) p(1)*exp(-((x-p(2)).^2)/(2*p(3).^2))-y;
    po(2) = find(y == max(y));
    po(1) = max(y);
    po(3) = 2;
    lb(1) = po(1)*.5;
    lb(2) = po(2)-5;
    lb(3) = .5;
    ub(1) = po(1)*2;
    ub(2) = po(2)+5;
    ub(3) = 4;
    f = lsqnonlin(F,po,lb,ub,opts);
    ep(n) = f(2);
end
return;


function ep = roughEdgeRefine(im,edgeEstimate,direction)
numframes = size(im,2);
ep = edgeEstimate;
buff = [2,15];
for n = 1:numframes

    switch direction
        case 1
            cur = im((edgeEstimate(n)-buff(1)):end,n);
        case -1
            cur = im((edgeEstimate(n)+buff(1)):-1:1,n);
    end
    %normalize linescan
    cur = cur-min(cur);
    
    thresh = max(cur)/2;
    ind = find(cur<thresh,1);
    if isempty(ind)
        continue;
    end
    switch direction
        case 1
            ep(n) = edgeEstimate(n) - buff(1) + ind - 1;
        case -1
            ep(n) = edgeEstimate(n) + buff(1) - ind + 1;
    end
end
    
return;




